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LATERAL  VELOCITY  VARIATIONS  AND  WATER 
CONTENT  ESTIMATED  FROM  MULTI-OFFSET  GROUND 

PENETRATING  RADAR 


Introduction 

Ground  penetrating  radar  (GPR)  is  the  recording  of  high  frequency  electromagnetic  waves 
that  have  reflected  from  subsurface  contrasts  in  dielectric  constant.  In  low-loss  media,  ground 
penetrating  radar  is  capable  of  producing  high  resolution  images  of  the  shallow  subsurface, 
compared  to  other  surface  electric  or  seismic  geophysical  techniques.  Most  commonly,  GPR 
data  is  collected  with  a  single  source-to-receiver  offset,  usually  a  minimum  offset.  In  the 
early  development  of  GPR,  multi-offset  common  midpoint  (CMP)  sounding  was  borrowed 
from  reflection  seismology  as  an  effective  technique  for  determining  glacial  ice  velocities 
(Gudmandsen,  1971).  As  the  use  of  GPR  expanded  to  rock  and  soil  surveys,  multi-offset 
radar  sounding  has  continued  to  be  used  primarily  for  velocity  sounding  at  one  or  just  a  few 
points  along  a  survey  line.  The  advantages  of  determining  velocity  with  this  method  are 
that  it  requires  no  prior  knowledge  of  the  subsurface,  is  not  intrusive,  uses  the  radar  data 
acquisition  system  only  and  can  determine  the  velocity  anywhere  within  the  survey.  The 
disadvantage  is  that  acquiring  multi-offset  data  with  current  GPR  systems  is  slow  (compared 
to  zero-offset  surveys),  or  impossible,  for  systems  with  fixed  source-receiver  offset. 

Recent  case  studies  have  shown  that  when  an  entire  GPR  survey  is  acquired  with  the 
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CMP  geometry,  multi-trace  reflection  seismic  processing  techniques  can  be  used  to  improve 
the  subsurface  radar  images  (Fisher  et  al.,  1992;  Gerlitz  et  al,  1993).  When  the  data 
are  collected  in  this  configuration,  normal  moveout  velocity  analysis  is  used  to  derive  a 
continuous  2-D  radar  stacking  velocity  field.  In  general,  lateral  radar  velocity  variation  has 
been  considered  to  be  small  within  the  spatial  range  of  most  surveys.  In  this  study,  we  show 
that,  in  fact,  lateral  variation  in  radar  velocity  can  be  significant  and  further,  that  with 
increasing  lateral  velocity  description  comes  improvement  in  the  results  of  multi-offset  radar 
processing. 

In  general,  ground  penetrating  radar  velocity  decreases  rapidly  with  depth  (Davis  and 
Annan,  1989).  This  is  primarily  a  result  of  increasing  water  content.  Topp  et  al.  (1980),  de¬ 
rived  an  empirical  relationship  between  radar  propagation  velocity  and  water  content.  Using 
this  relationship  we  can  estimate  water  content  from  the  interval  velocities  calculated  from 
normal  moveout  velocities.  This  interpretation  increases  the  potential  uses  of  radar  profiling 
in  ground  water  studies  and  contaminant  spill  monitoring.  Radar  interval  velocity  can  be 
used  to  estimate  water  content  when  the  subsurface  is  sufficiently  resistive  to  be  treated  as  a 
low-loss  media.  This  is  a  reasonable  assumption  where  radar  signals  penetrate  the  subsurface 
to  depths  of  meters  or  tens  of  meters.  In  partially  saturated  soils,  water  content  may  be 
interpreted  as  an  indicator  of  saturation.  In  fully  saturated  soils,  variations  in  water  content 
can  be  interpreted  as  variations  in  water  filled  porosity.  The  improvement  to  the  continuity 
and  depth-of-penetration  of  the  radar  image  combined  with  such  interpretations  of  the  ve- 
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locity  should  encourage  the  development  of  GPR  systems  capable  of  acquiring  multi-offset 
common  midpoint  data  efficiently  for  standard  surveying. 

Multi-offset  Data 

We  obtained  a  multi-offset  GPR  survey  from  Sensors  &  Software,  Inc.,  that  was  acquired 
at  the  Chalk  River  research  area,  operated  by  Atomic  Energy  of  Canada,  Limited.  The 
data  were  acquired  in  a  cooperation  of  Sensors  &  Software  Inc.,  the  Atomic  Energy  of 
Canada,  Ltd.  and  the  UT-Dallas  Geophysical  Consortium.  The  data  were  collected  with 
multiple  source  antenna  to  receiver  antenna  offsets,  such  that  after  re-arrangement,  1800 
CMP’s  spaced  every  0.25  m  were  defined.  The  data  were  acquired  using  a  pulseEKKO  IV 
digital  radar  system  with  100  MHz  antennas.  A  detailed  description  of  acquisition  geometry 
and  data  recording  parameters  can  be  found  in  Fisher  et  al.  (1992).  For  these  data,  each 
CMP  gather  has,  on  average,  ten  traces  with  offsets  ranging  from  0.5  m  to  20.0  m.  The 
standard  GPR  data  set  would  have  only  a  single  offset  trace  at  each  survey  station,  usually 
a  minimum  offset.  In  Figure  1,  the  minimum  offset  trace  profile  from  the  Chalk  River  data 
set  is  displayed.  This  shows  the  data  as  it  would  be  collected  in  a  single  offset  GPR  survey 
at  this  site.  The  subsurface  at  the  site  is  described  as  bedrock  covered  by  glacial  till  and 
fluvial  sand  deposits  (Davis  and  Annan,  1989;  Fisher  et  al 1992).  In  this  profile  we  see 
a  strong  channel  shaped  reflector  that  appears  to  be  the  depth-of-penetration  limit.  Also, 
there  is  a  clear  decrease  in  reflection  continuity  between  about  100-400  nanoseconds  from 
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NANOSEC 


Figure  1:  Minimum  offset  profile  extracted  from  multi-offset  GPR  data  at  Chalk  River. 
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CMP  900  to  CMP  1800.  Improvements  to  the  radar  image  in  these  areas,  in  particular,  are 

observed  in  the  final  CMP  processed  data  profiles. 

Figure  2a  shows  a  number  of  CMP  gathers  with  traces  arranged  in  offset  order  from  0.5  m 
to  20.0  m  within  each  CMP.  Some  of  the  approximately  hyperbolic  arrivals  from  reflections 
are  clearly  distinguished,  but  the  signal-to-noise  is  low  in  this  raw  data.  Pre-processing 
steps  were  applied  to  the  data  to  prepare  it  for  velocity  analysis  and  stacking.  The  data 
were  bandpass  filtered,  a  top  mute  (in  shot  gather  mode)  was  applied  to  remove  the  direct 
arrivals,  and  time  dependent  scaling  was  used  to  partially  correct  for  geometric  spreading 
loss.  A  filter  in  the  frequency  domain  was  applied  to  remove  some  spatial  aliasing  effects. 
AGC  (automatic  gain  control)  scaling  is  used  for  display.  Figure  2b  shows  the  same  CMP 
gathers  after  these  processes  have  been  applied.  These  filtered  data  are  used  in  the  normal 
moveout  velocity  analysis. 

Normal  Moveout  Velocity  Estimation 

Since  the  data  in  this  survey  were  all  collected  with  the  CMP  geometry,  normal  moveout 
(NMO)  velocity  analysis  can  be  applied  at  any  or  all  of  the  CMP’s  to  define  the  subsurface 
velocity.  How  many  velocity  analyses  are  required  depends  on  how  strongly  the  radar  velocity 
varies  laterally.  There  is  no  rule  or  formula  for  determining  an  optimum  number  of  velocity 
analyses  for  processing  a  CMP  data  set,  so  this  parameter  must  be  established  by  testing 
the  data  itself.  In  this  example,  we  show  how  increasing  the  number  of  velocity  analyses 
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Figure  2a:  Radar  CMP  gathers  in  offset  order  before  pre-processing.  Offsets  vary  from  0.5 
m  to  20.0  m  within  each  CMP. 
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Figure  2b:  Radar  CMP  gathers  after  pre-processing.  Offsets  vary  from  0.5  m  to  20.0  m 
within  each  CMP. 
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affects  the  final  CMP  stacked  image. 

There  are  a  variety  of  schemes  used  in  NMO  velocity  analysis  (Yilmaz,  1987).  We  have 
used  a  semblance  amplitude  approach  where  the  data  in  the  CMP’s  are  normal  moveout 
corrected  and  stacked  using  a  range  of  trial  velocities.  The  amplitudes  versus  time  over 
the  whole  range  are  then  contoured  and  displayed  as  a  velocity  spectrum.  The  peaks  of 
the  amplitude  mapping  are  chosen  to  define  the  1-D  velocity  function  at  each  CMP  being 
analyzed.  After  a  number  of  velocity  functions  have  been  defined,  a  2-D  velocity  profile 
is  created  by  interpolation.  For  a  statistical  method,  like  semblance  mapping,  the  more 
traces  there  are  in  any  individual  CMP,  and  the  larger  the  offset  range  (within  the  small 
spread  approximation),  the  better  the  resolution  in  the  velocity  spectrum.  In  practice,  the 
vertical  resolution  of  NMO  velocity  analysis  has  limits  such  that  only  the  highest  amplitude 
velocity  spectra  peaks  are  picked  for  the  traveltime  versus  velocity  function.  Therefore,  the 
NMO  velocity  function  will  usually  have  less  velocity  layers  defined  than  there  are  reflections 
observed  in  the  data. 

When  the  number  of  traces  per  gather  is  small  and  somewhat  noisy,  we  expect  the 
velocity  spectra  to  have  a  relatively  poor  signal-to-noise  ratio.  For  the  Chalk  River  data  we 
reduce  the  noise  in  the  velocity  spectra  by  combining,  and  sorting  by  offset,  the  traces  from 
twenty  CMP’s  centered  around  each  CMP  chosen  for  velocity  analysis.  A  typical  velocity 
spectrum,  for  a  single  CMP  (no  combination)  is  shown  in  Figure  3a.  When  the  spectrum  of 
the  combined  20  CMP’s  is  found,  as  shown  in  Figure  3b,  the  velocity  peaks  are  much  better 
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Figure  3:  (a)  Velocity  spectrum  of  CMP  755.  (b)  Velocity  spectrum  of  the  combination  of 
20  CMP’s  centered  on  CMP  755. 
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resolved  and  NMO  velocity  can  be  picked  with  more  confidence.  To  an  extent,  resolution 
in  the  spectra  seems  to  improve  with  increasing  traveltime.  This  is  due  in  part  to  the  fact 
that  for  near  surface  reflections,  the  further  offset  traces  are  muted  and  do  not  contribute 
to  the  analysis,  so  there  are  less  data  in  the  statistical  analysis.  Secondly,  for  a  constant 
test  velocity  interval,  changes  in  stacking  occur  more  slowly  as  velocity  increases,  leading, 
in  radar  data,  to  less  well  resolved  peaks  in  the  spectra  of  near  surface  reflections.  Overall 
width  of  the  velocity  spectra  peaks  in  Figure  3b  indicate  that  resolution  of  the  stacking 
velocity  is  on  the  order  of  +/  —  0.01  m/ns. 

Velocity  analyses  were  performed  at  regular  intervals  along  the  profile.  Figure  4  shows  the 
interpolated  NMO  velocity  fields  as  the  number  of  velocity  analyses  is  increased.  Figure  4a 
is  the  simple  flat  layered  profile  that  results  when  only  one  CMP  is  analyzed.  Velocity  varies 
from  about  .125  m/ns  near  the  surface  to  about  .06  m/ns  for  deep  reflectors.  The  deepest 
reflectors  for  which  velocity  is  estimated  is  at  about  750  ns  in  the  central  portion  of  the  data. 
Figure  4b  shows  the  2-D  velocity  profile  when  velocity  analyses  at  four  CMP’s  are  included. 
We  see  that  there  is  significant  lateral  variation  which  approximately  mirrors  the  sediment 
wedge  in  the  reflection  data.  Figure  4c  includes  velocity  analyses  at  nine  CMP’s,  Figure  4d 
at  eighteen  CMP’s  and  Figure  4e  at  35  CMP’s.  As  the  number  of  velocity  analyses  increase 
the  velocity  field  becomes  more  complicated  and  suggests  that  some  subsurface  properties 
not  directly  related  to  layering  are  affecting  the  velocity  structure. 
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Figure  4:  2-D  normal  moveout  velocity  profile  in  traveltime  after  analysis  at  a)  1  CMP,  b) 


4  CMP’s,  c)  9  CMP’s,  d)  18  CMP’s,  and  e)  35  CMP’s. 
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Stacked  Radar  Profiles 


Normal  moveout  corrections  calculated  from  the  NMO  velocity  field  remove  the  offset  de¬ 
pendence  of  the  reflection  traveltimes  such  that  the  data  can  be  treated  as  zero  offset  traces. 
After  the  correction  is  made,  the  traces  within  each  CMP  are  stacked  to  produce  a  single 
CMP  trace.  If  the  NMO  correction  is  done  with  the  correct  velocity  function,  the  stacked 
traces  have  an  improved  signal-to-noise  ratio. 

The  NMO  correction  is  based  on  the  assumption  that  the  subsurface  sampled  by  each 
CMP  can  be  adequately  modeled  as  a  sequence  of  horizontal  layers  with  uniform  interval 
velocity.  Steeply  dipping  reflectors  and  strong  velocity  gradients  test  the  applicability  of 
the  normal  moveout  and  stack  technique,  but  in  general  the  method  is  sufficiently  robust  to 
improve  data  quality  under  most  conditions.  If  lateral  velocity  variation  is  small  enough,  a 
single  NMO  velocity  function  can  be  used  to  calculate  the  moveout  correction  at  all  CMP’s. 
However,  where  lateral  velocity  variation  is  significant,  a  variable  NMO  velocity  field  defined 
by  functions  at  a  number  of  CMP’s  must  be  applied  to  obtain  the  best  results.  In  principle, 
an  NMO  velocity  defined  individually  by  velocity  analysis  for  every  CMP  should  yield  the 
most  accurate  result.  However,  a  spatial  limit  to  lateral  variation  in  NMO  velocity  is  usually 
reached  at  some  multiple  CMP  spacing. 

NMO  corrections  were  computed  for  the  Chalk  River  data  using  each  of  the  velocity 
fields  in  Figure  4.  In  Figure  5a  we  show  the  stacked  profile  when  NMO  is  defined  by  a  single 
velocity  analysis  (Figure  4a).  Figure  5b  is  the  stacked  profile  result  when  all  thirty-five 
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CMP  202  402  602  802  1002  1202  1402  1602 


Figure  5a:  CMP  stacked  radar  reflection  profile  using  velocity  profile  from  1  CMP,  corre¬ 
sponding  to  Figure  4a. 


13 


NANOSEC 


CMP  202  402  602  802  1002  1202  1402  1602 
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Figure  5b:  CMP  stacked  radar  reflection  profile  using  velocity  profile  from  35  CMP’s,  cor¬ 
responding  to  Figure  4e. 
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velocity  analyses  (Figure  4e)  are  used.  First  we  compare  these  to  the  standard  (no  stack) 
profile  shown  in  Figure  1.  The  stacked  data  profiles  show  a  strong  reduction  in  background 
noise  as  compared  to  single  offset.  Also,  some  deeper  reflections,  below  the  apparent  channel 
bottom  reflection,  have  been  greatly  enhanced,  in  particular,  the  reflector  at  about  600  ns 
from  CMP  250  to  CMP  1000.  Note  that  when  a  velocity  field  is  applied  the  exact  traveltime 
to  reflectors  changes  slightly.  This  is  a  consequence  of  the  normal  moveout  correction  which 
is  shifting  arrival  times  as  a  function  of  the  particular  velocity  field.  Comparing  the  single 
velocity  stack,  Figure  5a,  to  the  multiple  velocity  analyses  stack,  Figure  5b,  we  observe 
that  some  deeper  reflections  occur  down  to  750  ns  but  the  more  clear  improvement  is  in 
continuity  of  reflectors  throughout  the  section.  Overall,  stacking  with  even  just  one  velocity 
control  point  provides  the  major  increase  in  the  depth-of-penetration  over  no  stack.  With 
increasing  detail  in  the  velocity  field  there  is  some  increase  in  the  depth-of-penetration  but 
primarily  improves  reflector  continuity  and  time  structure.  The  fact  that  the  majority  of 
depth-of-penetration  increase  occurs  with  the  stack  from  even  a  single  velocity  analyses  as 
compared  to  no  stack,  suggests  that  it  is  the  multi-offset  nature  of  the  stack,  in  particular 
the  far  offsets,  that  provide  the  signal  from  the  deeper  reflectors. 

As  previously  noted,  the  data  between  about  200  to  400  ns  to  the  right  of  CMP  900 
shows  very  few  coherent  reflectors  in  the  standard  GPR  survey  shown  in  Figure  1.  In  fact 
the  left  half  of  the  profile  is  quite  different  than  the  right  half  of  the  profile.  With  the  NMO 
correction  and  stack,  the  region  on  the  right  changes  dramatically,  with  many  reflectors 
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emerging  from  the  stack.  When  only  one  velocity  is  applied,  Figure  5a,  this  area  of  the  data 
is  still  poorly  resolved.  However,  when  all  velocity  functions  are  included  in  the  velocity 
description,  the  stack  has  improved  continuity  of  a  number  of  reflections  from  the  left  side 
of  the  data  into  the  right  side.  Also,  note  that  the  deepest  reflection  on  the  right  side,  now 
appears  to  consist  of  a  series  of  step-like  events  which  may  be  interpreted  as  due  to  stream 
erosional  features.  Although  not  shown  here,  only  small  differences  between  the  stack  based 
on  eighteen  velocity  analyses  and  thirty-six  were  observed.  This  implies  that  lateral  velocity 
resolution  for  these  data  is  between  50  and  100  CMP’s  (12.5  to  25  meters). 


Radar  Propagation  Approximations 

Radar  propagation  is  fundamentally  limited  by  the  conductivity  of  the  subsurface  medium. 
Only  in  low-loss  media  can  radar  signals  penetrate  deep  enough  to  provide  a  useful  subsurface 
image.  For  radar,  low- loss  media  has  been  described  as  soil  or  rock  with  conductivity  less  than 
100  mS/m  (Davis  and  Annan,  1989).  Useful  approximations  for  describing  the  propagation 
of  radar  signals  can  be  found  by  considering  the  time  harmonic  solution  of  the  equation 


6E 

V2E  =  \mt  ——  -f  fi£ 
St 


S2E 
St 2 


(1) 


derived  from  Maxwell’s  equation  for  the  case  of  a  homogeneous  isotropic  medium.  E  is 
the  electric  field  intensity  and  the  constants  of  proportionality,  e,  /z  and  a  are  the  electric 
permittivity,  magnetic  permeability  and  conductivity  of  the  medium.  The  first  term  on  the 
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right  side  of  the  wave  Eq  (1)  represents  conduction  of  charge  due  to  the  applied  electric  field. 
The  second  term  describes  the  displacement  of  charge  due  to  the  field. 

A  solution  to  Eq  (1)  of  the  form 

E{z,t)  =  E0e~^zeiut  (2) 

yields  the  dispersion  relation 

k?  =  ificru  +  fieuj2.  (3) 

The  wavenumber,  &,  is  complex  and  can  be  written  as 

k  =  a  +  i/3  (4) 

where  the  attenuation  constant 

a  =  fl/f  (a/i  +  WS-i)  (5) 

and  phase  constant 

P  =  (\/l  +  ‘“2"«  +  J)  W 

are  both  real.  The  parameters,  c  and  Ke  are  the  electromagnetic  velocity  in  free  space  and 
the  real  part  of  the  dielectric  constant  of  the  medium.  Relative  susceptibility,  has  been 
eliminated  since  it  is  considered  to  be  unity  for  most  near  surface  earth  materials  (Telford 
et  al.,  1976).  The  solution  to  the  equation  can  then  be  written  as 

E{z,t)  =  E0e-aze~i^-pz)  (7) 
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which  is  the  expression  for  a  damped  plane  wave  propagating  with  phase  velocity 


(8) 


For  a  plane  wave  propagating  with  the  angular  frequency  u>,  the  ratio  of  conduction 
current  density  to  displacement  current  density  is  the  loss  tangent 


tan  8  =  — .  (9) 

ue 

In  materials  that  are  good  conductors,  displacement  currents  are  negligible  compared 
to  conduction  currents,  and  Eq  (1)  reduces  to  a  diffusion  equation,  i.e.,  the  fields  do  not 
propagate  as  electromagnetic  waves.  For  materials  with  low  conductivity,  and  when  the 
frequency  of  the  oscillating  electric  field  is  high  enough,  the  displacement  current  dominates 
over  the  conduction  current  and  electromagnetic  waves  will  propagate  (Stratton,  1941).  For 
GPR  systems,  which  by  definition  are  high  frequency,  the  loss  tangent  is  very  small  and  the 
diffusion  term  can  be  neglected.  In  this  case,  where  tan  8  <C  1,  the  phase  constant  reduces 
to 

P  ~  —  (10) 

c 


such  that 


and  or  reduces  such  that  a  depth  of  penetration,  dp ,  can  be  approximated  by 


(11) 


dp 


-  *  5  x  10-3  ^ 
a  a 


meters 


(12) 
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Since  common  soil  mixtures  have  dielectric  constants  less  than  fresh  water  («e  —  80),  only 
very  low  conductivity  materials  (<t  <  100  nriS/m )  will  allow  propagation  to  useful  depths. 
For  example  if  Ke  =  9  then  a  depth  of  penetration  of  one  meter  requires  that  a  =  15  mS/m. 

This  discussion  has  centered  on  the  time  harmonic  solution  to  the  wave  equation  and 
as  such  has  not  considered  that  the  material  properties  are  also  frequency  dependent  and 
therefore  should  be  treated  as  complex  numbers.  However,  Davis  and  Annan  (1989)  and 
others  have  shown  that  in  the  frequency  range  of  ground  penetrating  radar,  10-1000  MHz, 
low-loss  media  are  essentially  non-dispersive.  In  our  interpretation  we  treat  the  dielectric 
constant  as  real  and  related  to  the  propagation  velocity  by  the  simple  expression  found  in 
Eq  (11).  In  making  this  approximation  we  recognize  that  frequency  dependence  as  well  as 
many  other  factors  such  as  scattering  loss,  source/receiver  antenna  power  and  transmission 
characteristics,  and  ground  coupling  (Daniels,  1989)  are  not  being  accounted  for.  However, 
our  assumption  is  that  these  factors  change  much  more  slowly  than  dielectric  constant  if 
reflections  axe  observed  in  the  data. 

Interval  velocity  and  water  content 

To  interpret  the  NMO  velocity  field  derived  from  the  multi-offset  data,  it  is  necessary  to 
calculate  interval  velocities  and  find  the  relationship  of  radar  propagation  velocity  to  other 
geoelectric  properties.  Interval  velocity  utj„  was  calculated  using  the  Dix  formula  (Dix,  1955). 
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This  calculation  was  made  difficult  by  the  fact  that  radar  velocity  decreases  with  increas¬ 
ing  traveltime.  The  Dix  formulation  does  not  preclude  the  case  of  velocity  decreasing  with 
traveltime  but  we  found  that  with  decreasing  velocity  the  numerator  inside  the  square  root 
of  (13)  can  be  negative  if  the  traveltime  interval  is  small  or  the  NMO  velocity  change  is 
large.  Where  this  situation  was  encountered  we  used  the  average  velocity  of  laterally  adja¬ 
cent  intervals.  When  the  Dix  formula  is  applied  to  calculated  exact  NMO  traveltimes  and 
moveout  velocities  for  a  radar  velocity  model,  this  problem  is  not  observed.  This  implies 
that  the  Dix  inversion  when  applied  to  radar  data  is  very  sensitive  to  noise  in  the  velocity 
analysis.  In  consideration  of  this,  all  properties  calculated  from  the  interval  velocities  were 
subjected  to  strong  smoothing.  For  example,  time-to-depth  conversion  for  the  profile  was 
taken  to  be  the  least  squares  linear  fit  to  the  data  from  all  of  the  velocity  analyses. 

The  final  step  in  our  interpretation  scheme  is  to  relate  the  calculated  radar  interval 
velocities  to  water  content.  Interval  velocity  derived  from  reflection  moveout  curves  is  the 
group  velocity,  but  in  non-dispersive  media,  group  velocity  is  equal  to  phase  velocity  and  we 
can  use  the  approximation  (11)  to  convert  velocity  to  dielectric  constant.  Topp  et  al.  (1981) 
derived  an  empirical  relationship  between  measured  water  content  in  laboratory  samples  and 
dielectric  constant.  Their  empirical  equation  for  estimating  volumetric  wrater  content  0V  is 

ev  =  -5.3  x  10"2  +  2.92  x  10_2Ke  -  5.5  x  10_4/c2  +  4.3  x  10-6^  (14) 
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This  equation  is  most  appropriate  to  GPR  since  it  is  based  on  measurements  on  soil 
samples  of  varying  water  content  in  the  frequency  range  of  ground  penetrating  radar  systems 
with  a  laboratory  setup  equivalent  to  a  pulsed  radar  system. 

We  also  considered  two  other  functions  to  estimate  water  content.  The  CRIM  (complex 
refractive  index  method)  equation  is  used  for  interpretation  of  electromagnetic  propagation 
logs.  It  is  also  an  empirically  derived  mixing  law  relating  dielectric  constant  to  water  filled 
porosity,  (f>, 

<f>CRIM  = 

where  k  is  the  real  part  of  the  dielectric  constant  of  the  water  (kw)  and  rock  matrix  (/cOT) 
mixture  (Schlumberger,  1991).  To  apply  this  equation  to  GPR  data  we  assume  that  tan  8  <C 
1  for  all  components  of  the  mixture  so  that  the  equation  reduces  to 


y/l— tan2  | 


-  y/i^ 


y/  *>w 


y/l— tan2  | 


\ZKm 


(15) 


4>crim 


y/H  —  y/K^ 


(16) 


y/  ■y/^'771 

Another  mixing  formulation  is  discussed  by  Sen  et  al.  (1981)  for  estimating  water  sat¬ 
urated  porosity.  This  is  the  two-phase  Hanai-Bruggeman  equation  which  in  the  frequency 
range  of  radar  can  be  approximated  by 


)i(JLr5=-)  (17) 

\  K  /  \  u  fcm  f 

where  m  in  the  exponent  is  the  geometric  factor  from  Archie’s  law  (Samstag,  1992).  This 
factor  is  most  commonly  taken  to  be  m  =  2  (Sheriff,  1991)  but  in  Jackson  et  al.  (1978)  it 
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is  empirically  shown  to  range  in  marine  sands  from  m  =  1.52  for  shaly  sand  to  m  =  1.9  for 
platey  shell  fragments.  In  saturated  soils  the  porosity  calculated  from  (16)  or  (17)  is  equal  to 
the  fractional  water  content  of  Eq  (14).  Figure  6  shows  a  comparison  of  the  the  water  content 
versus  dielectric  constant  calculated  with  these  equations.  For  equations  (16)  and  (17)  we 
used  kw  =  80  and  /cm  =  4.5  to  approximate  a  fresh  water  saturated  sand.  The  geometric 
factor  in  Eq  (17)  was  taken  to  be  m  =  2  for  simplicity  since  we  are  measuring  quantities 
on  a  gross  scale  compared  to  the  micro-geometry  reflected  in  the  m  factor.  From  the  plot 
we  can  see  that  any  of  these  relationships  would  yield  similar  results  in  converting  dielectric 
constant  to  water  content  but  the  point  of  this  calculation  is  to  show  that  the  Topp  Eq  (14) 
is  in  good  agreement  with  other  estimation  techniques.  If  the  GPR  data  was  collected  at 
multiple  frequencies  it  would  be  more  appropriate  to  use  the  Hanai-Bruggeman  equation 
which,  in  its  complete  form,  is  frequency  dependent.  It  should  also  be  noted  that  the  actual 
measured  property,  interval  velocity,  changes  more  slowly  with  water  content  as  velocity 
decreases.  Since  radar  velocity  in  general  decreases  with  depth,  decreasing  sensitivity  of  the 
inversion  to  variations  in  water  content  is  expected  to  occur  at  greater  depths. 

We  used  the  Topp  equation  to  make  the  estimate  of  water  content  shown  in  Figure  7 
from  the  interval  velocity.  In  order  to  compare  this  plot  to  the  reflection  profiles,  the  water 
content  is  plotted  linear  in  time,  but  with  the  approximate  depths  indicated  on  the  right. 
The  water  content  estimates  were  binned  into  20  ns  intervals  and  then  smoothed  three 
times  with  a  nine  point  averaging  window.  The  result  shows  clearly  a  zone  of  increasingly 
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Figure  6:  Relating  dielectric  constant  to  water  content  in  low-loss  media.  Comparison  of 
Topp  et  al.  (1981)  empirical  relation  to  the  CRIM  (Schlumberger,  1991)  and  Hanai- 
Bruggeman  (Sen  et  al.,  1981)  equations. 
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Figure  7:  Water  content  estimated  from  the  ground  penetrating  radar  interval  velocity  using 
the  Topp  et  al.  (1981)  equation.  Values  have  been  smoothed  over  150  CMP’s  laterally 
and  120  nanoseconds  in  time. 
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shallow  high  water  content  from  left  to  right  along  the  profile.  This  can  be  interpreted  as 
an  indication  of  a  rising  water  table.  Since  this  region  of  high  water  content  cuts  across  the 
detailed  reflection  structure,  it  implies  that  water  filled  porosity  and  permeability  pathways 
are  not  constrained  to  apparent  stratigraphic  structure.  Lateral  variations  in  water  content 
will  occur  within  depositional  units  as  sand  and  clay  ratios  and  grain  size  vary.  Since  water 
distributed  across  stratigraphic  units  is  not  likely  to  be  connate  water,  it  is  reasonable  to 
conclude  that  the  zones  of  high  water  content  are  also  the  zones  of  highest  permeability. 
There  is  also  some  indication  that  there  are  areas  of  decreasing  water  content  at  depth 
which  can  be  interpreted  to  be  a  result  of  differences  in  soil  porosity  determined  by  soil  type 
and  compaction. 

Conclusion 

Ground  penetrating  radar  surveys  collected  with  the  CMP  multi-offset  geometry  yield  im¬ 
proved  subsurface  images  over  single  offset  surveys.  We  have  shown  that  the  CMP  profile 
is  itself  improved  as  the  number  of  velocity  analyses  is  increased.  This  leads  to  the  con¬ 
clusion  that  lateral  variation  in  radar  propagation  velocity  can  be  significant  even  over  the 
limited  range  of  a  typical  GPR  survey.  The  CMP  stacking  process  yields  improved  depth- 
of-penetration  over  a  single  offset  survey  while  detailed  velocity  analysis  yields  improvement 
to  continuity  of  reflectors  throughout  the  stacked  profile. 

In  this  study  we  have  also  attempted  to  connect  the  practical  measurement  of  radar 
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velocity  from  the  multi-offset  data  with  the  theoretical  and  laboratory  relationships  between 
water  content  and  dielectric  constant.  We  have  shown  a  practical  approach  to  estimating 
water  content  from  these  velocities  based  on  the  assumption  that  media  conducive  to  radar 
propagation  is  essentially  non-dispersive.  The  results  of  applying  this  method  to  the  Chalk 
River  data  illustrates  the  method  but  due  to  the  lack  of  comparative  data,  in  particular  field 
measured  water  content  or  dielectric  constant,  we  cannot  establish  how  accurate  these  results 
are.  It  should  be  considered  for  such  approximations  that  the  accuracy  is  no  better  than  that 
of  the  velocity  which  is  estimated  to  be  within  10-20%  of  the  average.  Some  improvement 
to  this  could  be  made  by  increasing  the  number  of  source  to  receiver  offsets  recorded  for 
each  CMP%.  We  also  suggest  that  a  test  of  this  interpretation  method  should  be  done  by 
measuring  subsurface  water  content  in  an  area  while  acquiring  a  multi-offset  radar  survey. 
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